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Abstract 

For the symmetric case of space- fractional diffusion processes (whose 
basic analytic theory has been developed in 1952 by Feller via in- 
version of Riesz potential operators) we present three random walk 
models discrete in space and time. We show that for properly scaled 
transition to vanishing space and time steps these models converge in 
distribution to the corresponding time-parameterized stable probabil- 
ity distribution. Finally, we analyze in detail a model, discrete in time 
but continuous in space, recently proposed by Chechkin and Gonchar. 
REMARK: Concerning the inversion of the Riesz potential operator 
/q let us point out that its common hyper-singular integral repre- 
sentation fails for a = 1. In our Section 2 we have shown that the 
corresponding hyper-singular representation for the inverse operator 
Dq can be obtained also in the critical (often excluded) case a = 1, 
by analytic continuation. 
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1 Introduction: concepts and notations 

By a " space- fractional" diffusion process (or Levy- Feller diffusion process) 
we mean a process of diffusion of an extensive quantity with density u{x,t) 
governed by an evolution equation 

(1.1) = D^u{x, t) {xeR,t> 0) 
with an initial condition 

(1.2) uix,0) = f{x) (/GLi(R)). 

Wc interpret x as space, t as time variable. Dg is a pseudo-differential 
operator acting with respect to the space variable x, its symbol being 

(1.3) Df{K) = -|«;|«e^(^ig^'^)^^/2 = -|«;|«i(sign«;)^^ 
and the real parameters a and 9 are restricted by the inequalities 

(1-4) < a < 2, \0\< { ^ .r ~ n 

^ ^ ''-l2-Q;ifl<Q;<2. 

For a sufficiently well-behaved function or generalized function defined on 
R we denote by its Fourier transform: 

+00 

(1.5) $(k) = J e^*^^ 0(x) dx (k e R) . 

—oo 

Then, for a generic linear pseudo-differential operator A acting on these func- 
tions, its symbol A{k,) turns out to be defined through the Fourier represen- 
tation of {A4>) (x) , namely A<j){K) = A{k) 4>{k) . An often applicable practical 
rule is 

(1.6) A{k) = (^e-*'^^)e^'^^ , (/« e R) . 

If B is another pseudo-differential operator, then we have A B{k) = A(k,) B(k) . 
Let us remark, that we chose (1.5) to define the Fourier transform in agree- 
ment with the common terminology of probability theory. 

Introducing the stable probability density Pa{x; 6) whose characteristic func- 
tion (Fourier transform) is 

(1.7) P^{k; e) = exp ( - |/^re^(«ig^ '^)^^/2) {k e R) 
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and rescaling Pa{x;9) for a; e R, t > by the similarity variable xt ^ we 
obtain the time-dependent stable probability density 

(1.8) ga{x,t;e) =t~^Pa{xt-^;9) (x G R, t > 0) 
with which we can write the solution to (1.1), (1-2) in the form 

+ 00 

(1.9) uix,t)^ I 9aix-U;0)fiOdC- 

—oo 

Then, for alH > 0, we have 

(1.10) < 



+ 00 



u{-, i) e C°° n Li(R) , / u{x, t)dx^ J f{x) dx 

— oo — oo 

, f{x) > for all X e R =^ u{x,t) > for all x e R. 



If the initial function is a probability density then so is also the function 
u{x,t) and we have in (1.1), (1.2) the description of a Markov process. 

For orientation on the general theory of stable probability distributions we 
recommend in particular [F71] and [F52], Feller's parameterization being 
close to ours. In our foregoing considerations we essentially have surveyed 
results of Feller's pioneering paper [F52]. For a few parameters pairs (a, 6) 
representations of ^^(a;; 6) in terms of elementary or well- investigated special 
functions are available in the literature but in other parameterization (see 
e.g. [Zo], [Sc], [SaT]). We here content ourselves with recognizing the classical 
Gauss process and the Cauchy process, respectively, in 

1 1 t 

,,(x, t; 0) = ^t-l exp ( - -) , ,,(x, t; 0) = --^ . 



Feller in [F52] has shown that Dq (in case a ^ 1) can be viewed as inverse 
to the (later in [SKM] so called) Feller potential operator which is a linear 
combination of two Weyl integral operators. In the sequel we will restrict 
attention to the symmetric case ^ = retaining however, the index in 
order to be in concordance with the notation of our previous papers [GM98], 
[GFM], [GM99]. In [GM98] we have discussed for a 7^ 1 a random walk 
model for the whole range of values 9] here we now concentrate on models 
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not treated there, requiring however ^ = . This means that henceforth we 
will treat the evolution equation 



(1.11) ^ ^^^^ {xeR,t> 0) 

with an initial condition 



(1.12) u{x,0) = f{x) (/eLi(R)). 
The symbol of the pseudo-differential operator Dq is 

(1.13) d^{k) = -\Kr, 

and the fundamental solution to (1.11) (namely, for u{x, 0) = 5{x) — Dirac's 
delta function) is the function u{x,t) — ga{x,t;0) whose Fourier transform 
is 



(1.14) ^(«;,t;0) = exp( -i|«;|°^) («;eR,t>0). 
Then, see (1.7) and (1.8), we get by the Fourier inversion formula 

(1.15) ga{x,t;0)^^ /^"^ e"*^^ exp(-t|fi;|") = r * p«(a;t-- ; 0) . 

ZTT J-(X} 

We will denote by S{t) for i > the random variable whose probabihty 
density is given by Qaix, t; 0) . 

Remcirk 1.1. We call a process described by (1.1), (1.2) a "space fractional" 
diffusion process because it is a "fractional" generalization of the classical 
diffusion process which is recovered by taking a = 2 , 6 = . Writing (1.13) 
in the form Dq{k) = — (k^)"/^ and observing that the operator defined 
by (D2 0) (x) = , has the symbol -k^ , we see that = - {-D^f'^ , 
hence the operator is the negative of a fractional power of the (positive 
definite) operator —D"^. By calling a process described by (1.1), (1.2) also 
a Levy-Feller diffusion process we honour both Levy and Feller for their 
essential contributions [L25], [L54], [F52]. 

Our aim in the following sections is to derive discrete-space discrete-time 
random walk models approximating the space-fractional diffusion process 
(henceforth considered as a Markov process) described by (1.11) with u{x, 0) — 
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5{x) . We shall show that for properly scaled transition to vanishing space and 
time steps there is convergence in distribution to the probability distribution 
whose density is ga{x,t;0). We shall give heuristic motivations for choos- 
ing our concrete models, using in a formal way calculations with symbols. 
The lack of rigour in these derivations will hopefully not be too annoying to 
the pure analyst, it will be remedied in the final end by rigorous proofs of 
convergence. 



2 Operators and symbols 



In this section we give a survey on the relevant operators and their symbols 
thereby always assuming 



(2.1) 



< a < 2. 



For general orientation and more rigorous treatment we refer to [SKM] , [R] 
and [F52]. We need the operator D of differentiation, the Weyl operators 
and their (formal) inverses /^", /Z", the Riesz potential operator 
/q whose negative inverse — /(f" is our pseudo-differential operator Dq (if 
q; 7^ 1), and the Hilbert transform operator H. For sufficiently well behaved 
functions defined on R and if required with appropriate understanding of 
the occurring integrals as Cauchy principal values we have (for x e R) 



(2.2) 



(D</.)(x) = ^0(a:) = 0'(x), 



(2.3) 



(/«(/>) (x) 



1 



1 



(2.4) 
(2.5) 



±DI]r'^ if < a < 1 , 

^2j2-a jf 1< a < 2 , 



(7^ <f>){x) + {i^ 0)(x) \x - er-^ m 



2cos(Q;7r/2) 



2r(Q;) cos (Q;7r/2) 



for a 7^ 1, 
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(2.6) ^H<P){x) = - r°°^d^. 



Recalling formula (1.6) for calculation of symbols we get by direct computa- 
tion 

(2.7) D^-iK, 



(2.8) H^isiguK. 

We take the symbols of the Weyl operators and the Riesz potential operators 
from [R, Theorem 4.10] as 

(2.9) = (tz/^)-" = \^\-a^±t{signK)a7T/2 ^ |^|-a ^ ±a sign ^ 
from which by addition we get 

(2.10) = . 
As already remarked and as used in [F52] 

(2.11) = -/q-" for a^l 
in agreement with the property of symbols 

(2.12) D^{k) = -IkI'' = -(J^iK))'^ for a 7^1. 
In the special case a = 1 we observe that by (2.7) and (2.8) 

(2.13) dI{k) = -\k\ = -D{k)H{k), 

hence as already observed in [F52] (with a sign-modified version of the Hilbert 
transform) 

(2.14) Dl^-DH. 

Let us yet exhibit another representation of our pseudo-differential operator 
Dq . Via the semi-group property (see [F52] or [SKM]) 

I^I^^lt^ if «,/3e(0,l), « + /3<l, 



analytic continuation to negative exponents can be justified and thus from 
(2.5) and (2.11) tlie formula 



In Section 4 we shall find random walk schemes in the case a 7^ 1 by approx- 
imating in (2.15) the operators and /Z" with the Griinwald-Letnikov 
discretization. In Section 5 a second random walk scheme, for the whole 
range < « < 2, will be obtained by a straightforward discrete approxima- 
tion of hyper singular integrals for Iq"^ and DH. Prom [SKM, formula (12.1')] 
we take (for 0<q;<2, a^l) 

(Albj (io '^^W-2r(_«)cos(a7r/2)yo 

Quite formally we can obtain (2.16) by replacing in (2.5) the integrand by 
- 0, then replacing a by -a, splitting = + j^, here 
regularizing the right hand side hypersingular integrals by subtracting (f>{x) 
in the numerators, finally substituting — ^ for ^ in the first right hand side 
integral and then putting both integrals together. 

For convenience we simplify the coefficient in (2.16), introducing 

(2.17) 6(a) := = - r(a + 1) sin (a7r/2) , 

^ ' ^ ' 2V{-a) cos(a7r/2) tt ^ ' ^ ' 

where for the latter equality we have used sin(Q;7r) = 2 sin {an jl) cos {an /2) 
and the reflection formula for the gamma function T{—a)V{a+l) = —n/ sin(Q;7r) 
Then, for < a < 2, a 7^ 1 , 
(2.18) 

Note that for (f) e C^(R) and (p{x) bounded the integral is finite as an 
improper Riemann integral and observe that b{a) > for the admitted values 

of a. 

Because (2.17) gives 6(1) = l/vr we are tempted to believe the formula 

(2.19) = i r ^ 

TT Jo 
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We will use (2.18) and (2.19) as motivation for our second random walk 
scheme in which the parameter value a = 1 does no longer play a special 
role. We can indeed obtain (2.19) by looking at (2.14), then formally differ- 
entiating (behind the integral sign) (2.6) and then splitting and regularizing 
the resulting hypersingular integral in the same way as we have done in the 
case a 7^ 1 . 

In view of (2.18) and (2.19) we now have, with b{a) given in (2.17), 
(2.20) 

(CWM^6(a)f + + for 0<a<2. 

Unfortunately, this formula loses its meaning in the case a — 2 . 

3 General structure of the random walk mod- 
els 

What do we mean by a random walk model, discrete in space and discrete in 
time, for a Markov process"? Let us be given a random variable Y taking its 
values in the set Z of integers according to the probabilities 

(3.1) P{Y^k)^pk (keZ) 
with 

(3.2) all Pk > and = 1 . 

keZ 

We discretize the space variable a; G R and the time variable t > by grid 
points Xj = jh and instants tn = nr, with /i > 0, r > 0, j e Z, n e Nq. 
Then, defining random variables 

(3.3) Sn^J2ihYm)^hJ2Ym (neN) 

m=l m=l 

with the Ym as independent identically distributed random variables, all hav- 
ing the same probability distribution as the random variable Y, we interpret 
Sn as the position at time t„ of a random walker starting in the point x — 
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at t — . Denoting by yj{tn) the probability of sojourn of the walker in point 
Xj at instant tn, the recursion S'„+i — Sn + hY^ implies 

(3.4) yjitn+i) = PkVj-kitn) U e Z, n e No), 

ikeZ 

and the walker starting at point xq = means yo{0) = 1 and yj{0) = for 
j ^ . However, in the recursion scheme (3.4) it is legitimate to use a more 
general initial sojourn probability distribution {yj{0)\j £ Z). 

There is yet another possible interpretation of (3.4), namely as a scheme 
of redistribution of an extensive quantity (e.g. mass, charge, in the random 
walk interpretation probability), yjitn) being considered as a clump of this 
extensive quantity sitting in point Xj at instant t^, . Then (3.4) describes a 
conservative and non-negativity preserving redistribution scheme. In fact, 
for all n e N it follows from (3.2) that, in analogy to (1.10), 

E%(^n) = E%(o) if EI%(o)l<°°' 

jeZ jgZ jeZ 

all yj{tn) > if all yj(0) > 0. 

Such discrete redistribution schemes have been used by one of the authors in 
discretization of diffusion processes governed by second order linear parabolic 
differential equations ([G70], [G78], [GN]) as they discretely imitate essential 
properties of the continuous process. 

We come nearer to the Cauchy problem (1.11), (1-12) by intending yj{tn) as 
approximation to 

Xj+h/2 
j u{x,tn)dx 

Xj-h/2 

which, if u{-,tn) is continuous, is ~ hu{xj,tn) for small h. 

We want to show that for proper choice of the probability distribution of the 
random variable Y and well-scaled transition 

(3.5) T = a{h) , (7 strictly monotonic, a{h) as h ^ 

the random walk "converges" in some sense to the Markov process with den- 
sity u{x,t) described by (1.11), (1-12) in the case that the initial function 
is a probability density. More specifically, we will prove for fixed t > 0, 
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t = tn = nr with N9n = t/r— >oo (and proper scaling of h and r) that 
the random variable Sn of (3.3) converges in distribution (other terminology: 
in law) to the random variable S{t) whose density is (7Q,(-,t;0) , the funda- 
mental solution (1.15) of (1.11). Observing that the distribution function 
Ga{x, t; 0) = daii, t] 0) is continuous in x (due to the fast decay in |«;| 
of ga{n,t;0) — exp(— i|«;|") the density ga{-,t;0) is in C°°(R)) and invoking 
the continuity theorem of probability theory (see, e.g., [B, Theorem 8.28]), 
all we have to do is to show that for all k G R the characteristic function 
y{K,, t; h) of the random variable Sn tends to exp(— t|fi;|") as h ^ . Note the 
equivalences following from t — tn — nr and (3.5) 

(3.6) n^oo<^/i^O<^r^O. 

The general form of the characteristic function y(«;, h) can be found via the 
generating functions 

(3.7) p{z) = J2 Pj^^^ y{z,tn) = J2 VjiQ 

As probabilities both the "Pj and yj{tn) sum up to 1 if added over the index j, 
hence these series converge absolutely and uniformly on the periphery 1 2; | = 1 
of the unit circle, and so the functions p and yn are there uniformly contin- 
uous. The random walk Sn starting at a; = 0, we have (using the Kronecker 
symbol) yj{0) — Sjo, and the recursion (3.4) being a discrete convolution we 
get 

(3.8) y{z) = {p{z)r. 

Replacing in (3.7) by e'^'^^^ = e*"-''' we obtain the corresponding character- 
istic functions («; e R) 

(3.9) p(k; h) = p{e'^^) , y(«:, h) = y{e'^h , tn) = (pie'^^J . 

Recalling our fixation of t = t„ = nr = na{h) > , the scaling relation (3.5) 
and the equivalences (3.6) we have to show that 

(3.10) y(«;, /i) — >■ exp(— i|«;|") for n — > 00 , 
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or, equivalently 

(3.11) ^— logp(e^''^) ^ -|/«|" as h^O. 

a{h) 

In the following sections we shall exhibit (3.10) as true for specific choices of 
the probabilities pj and scalings r = ct(/i) . The fact that, strictly speaking, 
r = t/n = (j{h) and h in (3.6) and (3.11) arc running through discrete sets 
will turn out as irrelevant for the proof of (3.11). 



4 The Griinwald-Letnikov random walk 



An idea suggesting itself is to discretize in (1.11) the time derivative ^ by 

a two-level difference quotient and the operators /jp" and Jz" (see (2.15)) 
by the Griinwald-Letnikov approximation (see, e.g. , [SKM], [P]). This idea 
leads to 

with the operators still to be specified. We must exclude the singular 

case a = 1 , hence will distinguish from now on the cases 

(a) < a < 1 , (b) l<a<2. 



We define 

(4.2) hl±''yj{tn) = h-^ f:(-l)" (^) Vj^kitn) in case (a) , 

(4.3) /./r = h-" f;(-l)' (^) Vm^kitn) in case (b) . 

Note in case (b) the shift of the index, that is required in order to obtain 
non-negative values for all transition probabilities pj . 
Solving (4.1) for yj{tn+i) , thereby scaling by 

(4.4) T = i^h"' =: a(h) , 
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gives (remember (3.4)) 
(4.5) 

with in case (a) 



feeZ 



(4.6)po = l Pk-{-lf\^'-^ for MO, 

cos(Q;7r/2) 2cos(Q;7r/2) \^|fc|y 



in case (b) 



(4.7) 



Po = 1 + 



P±i = 



cos (a7r/2) \1/ ' 
/i 



P±k = (-1) 



2 cos (a7r/2) 



a \ 












1 + 


(3 






a 



2cos(a7r/2) \k + 1 
Then, all pjt > if 

(4.8) < < cos(Q;7r/2) in case (a) , 



for A; = 2,3,... 



(4.9) 



< /J, < \ cos(Q;7r/2)|/Q; in case (b) . 



Note that cos(Q;7r/2) < in case (b). In both cases, by rearrangement of 
series. 



E(-i) 



^ cos (a7r/2) 

+00 



1-0 = 1 



Remcirk 4.1. For all a > the series '^{—lY ( j) absolutely convergent 

because = 0(j~°~^) for j — )■ cx) . 

We see that, under the conditions (4.8) or (4.9), respectively, we can put 

P(Y^k)^Pk (keZ) 

for the random variable Y of (3.1). Using (4.6) and (4.7) we identify the 
generating function p of (3.7) as 



(4.10) p{z) = 1 



2cos(Q;7r/2) 



{{l-z)" + {l-z-T} in case (a). 
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(4.11) p{z) - 1 - ^^^^ ^^^^^^ {z-\l - zT + z{l - z-'T} in case (b). 



Let us verify the limit relation (3.11) which implies (3.10). Because of the 
symmetry relation (for z — e*'^'^, k e R) 

p{z)^p{z~^) implying p(e"^'*) = p(e""^'*) 

it suffices to verify (3.11) for k > (the special case k — Q being trivial). 
Let be «; > 0. Then in case (a) we have 

^ cos(a7r/2) ^ ' ; v ; 

since (1 — zY ~ {—tKh)" = e"*"''/^ (Kh)" . With the scaling (4.4), namely 
a{h) = fih"", follows (3.11). 

In case (b) we have 

P{z) = 1 J^-— ^{z-\l - z^) , 

' cos(a7r/2) ^ ^ ' " 

and an analogous calculation gives 

1 1 ~t iKh\ c.cos[(Q;7r/2) + K/i] 

■ logp(e ) -k" — / \ — - as /i ^- , 



a{h) cos(Q;7r/2) 
hence again (3.11). As result wc have 

Theorem 4.2. Distinguish the cases (a) < a < 1 , (6) 1 < a < 2 . 
Define the probabilities Pk — P{Y — k) in case (a) by (4-6) with restriction 
(4-8), in case (b) by (4-V with restriction (4-9). Let the scaling relation 
T = fi = a{h) hold and let for fixed t > the index n = t/r run through N 
towards oo. Then the random variable Sn of (3.3) converges in distribution 
to the random variable S{t) whose probability density is given by (1-15) as 
ga{x, t; 0) . 

Remcirk 4.3. In the special case a — 2 the familiar explicit difference 
scheme 

Vjitn+l) = (1 - 2/x) yj{tr^) +liyj{tn-l) + //%(i„+i) 

is recovered from (4.7), and (4.9) goes over into the well-known stability 
condition < < 1/2 . 
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Remark 4.4. The case a = 1 is singular. For o; — > 1 both upper bounds in 
(4.8) and (4.9) tend to 0, and the denominators occurring in the definitions 
of the probabihties pk tend to zero. 

Remark 4.5. A motivation for the Griinwald-Letnikov approximation of 
can be drawn from the fact that z = e*"'^ is the symbol of the backward 
shift by a step h: With 

{Th(P) (x) = (P{x + h) , {T_h(t>) {x) = (t>{x - h) , 

we have 

Prom the symbol h~^{l — z)— hD+{K) of the usual backward approximation 

{hD+<P){x) = h-\ct>{x) - ct>{x - h)) = h-' Umx) 

we arrive by analogy at the symbol h~'^{l — z)°' as a candidate for the symbol 
of the operator . Analogously we get h~°'{l — z~^)°' as the symbol for 

the operator /j/Z" ■ We use the corresponding approximations in case (a) 
< a < 1 . 

In case (b) 1 < a < 2 we use the form /^^" of the Riemann-Liouville left 
inverse of the operator 7^, and put ^7^" = hI+~°' ■ The corresponding 
symbol then is, with symmetrically {hD'^(t>){x) — h"'^ (0(x + h) — 2(f){x) + 
Hx - h)) , 

„I^^{k) = hD\k) = (^"' - 2 + ;^) /i'-" (1 - z)-^^-'^^ 

= z-^ (1 - 2^ + z'') (1 - zy-'' = /i-" z-^ (1 - zY . 

The symbol h^'"' (1 — z)°'~'^ for /j/^"" here has been derived by the formal 
stipulation ^7^"" = /i-D+*-^~°^^ using hD+{K) — {1 — z) . 

Analogously we get h""' z {1 — z"^)" as symbol of ^71"^ . 

Remark 4.6. In [GM98] and [GM99] we have exploited the Griinwald- 
Letnikov random walks in the more general setting of not necessarily symmet- 
ric Levy- Feller diffusion (see Section 1). The proof for the case of symmetry 
{9 = 0) given in the present paper is considerably simpler. 
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5 The Gillis- Weiss random walk 



Gillis and Weiss in 1970 (see [GiW]) showed (we interpret one of their results 
in the language of probability theory) that every symmetric random variable 
Y with values in Z and asymptotically P{Y = A;) ~ c/\k\°''^^ (where c > 
0) lies in the domain of attraction of the corresponding symmetric Levy 
distribution, hence can be used for an approximating random walk in the 
sense of Section 3. Only assuming their asymptotics they naturally cannot 
describe precisely how the coefficients fj, and A of the scaling law appear in 
the transition probabilities. However, from their analysis we can deduce that 
the scahng law is of the form 

r = a{h) = nh" if 0<a<2, t = a{h) = Xh^\logh\ if a = 2 . 

Remarkably, the parameter value o; = 1 is not singular, but the scaling law 
becomes discontinuous at « = 2, thus giving an example of a distribution with 
non-finite variance lying in the domain of attraction of the normal (Gauss) 
distribution. 

We will now re- work and complement their analysis in the framework of our 
Section 3 for the special symmetric probability distribution {pk\k & Z) with 



The parameter a is only restricted as in (1.4) by < a < 2 . Differently from 
Gillis and Weiss we motivate this choice of probabilities by (2.20), where the 
special character of the value a — 2 already becomes visible. So, assume 
meanwhile < a < 2 . 



Discretizing D^u via a straightforward quadrature formula for the right hand 
side of (2.20) as 



oo 



(5.1) 



po = 1 - 2A 5] , pfe = A|A;|-("+^) for k , 



k=l 



where (so that po > 0) A is restricted by 



(5.2) 




(5.3) 



h 



D^yjitn) = b{a) h Y: 



k=l 



oo 



yj+k{tn) - 2yj{tn) + Uj-kitn) 
{kh)^+^ 
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and solving the equation 



yj{tn+l) - yjjtn) ^a^, x 

T 

for yj{tn+i) we identify the transition probabihties pk in (3.4) as 

(5.4) po^l-2fxb{a)C{a + l), Pk^ f^b{a)\k\-^'^+^^ for k , 
with /X = /i~"r , b{a) = T{a + 1) sin {aiT/2)/iT and the Riemann ('-function 

oo 

(5.5) C{z) = J2 for > 1. 

k=l 

Obviously Y^^eT^Pk ~ 1 ? ^'^^ the non-negativity condition in (3.2) requires 

1 TT 



(5.6) < < 



2 6(a) C(a + 1) 2T{a + 1) sin (a7r/2) C(a + 1) 



We want to free the parameter value a = 2 from its singular character. 
Recalling (2.17) we see that 6(2) = , so that in (5.4) po — 1 and all Pk — ^ 
for A; 7^ whereas the upper bound for in (5.6) tends to oo as a — > 2—. 
This degenerate random walk obtained in (5.4) by formally setting a = 2 
being neither interesting nor useful we replace fib{a) by A and obtain the 
transition probabilities in the form (5.1) with restriction (5.2). In (5.1) the 
special value a — 2 seems to be a quite regular value, and we shall see that 
we have a valid random walk model for all a obeying < a < 2 . However a 
price must be paid. Whereas for < a < 2 we can scale by r = /x /i" we can 
no longer do so in the case a = 2 . So, assume henceforth (if not explicitly 
stated otherwise) the condition (5.2). 

We have now the generating function 

oo 

(5.7) p{z) = 1 - 2XC{a + 1) + A ^ k-^"+'^ {z'' + z-'') 

k=l 

with z — e*'*'*, K e R. With the polylogarithmic function 
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we can write 



(5.8) p{z) = 1 - 2XC{a+ 1) + \ mz,a + 1) + <^{z-\a + 1)} 

and could carry out the required asymptotic analysis by specializing some of 
the formulas in [T] . See also [EHTF] and [Le] for properties of the polyloga- 
rithmic function and the more general Lerch function. We prefer, however, 
the direct way to obtain (3.11). This asymptotic relation is trivial for k = , 
and because of p{e^'^^) — p{e"^'^^) , it suffices to treat the case k > what we 
now will do. 

Prom the common integral representation of the gamma function we take 



1 r'^ 
rfa + l h 



V{a + l) Jo 
and using z"^ — z we get 

(5.9) p{z) ^ 1 - 2X^'y{z) 

with 
(5.10) 

/•CO °° rOD p—u 1 ~ 



The last equality in (5.10) has been obtained by summing the two involved 
geometric series. 

In the Appendix wc have performed in detail the required asymptotic analysis 
of 3?7(z) for u = K,h ^ 0+ (/t fixed), which is resumed in formulas (A. 6) and 
(A. 7). Insertion of these asymptotic behaviours into (5.9) yields 

(5.11) logp(e-'^) ^ - In^h" if < a < 2, « ^ , 

r{a + 1) sm (oiTr/z) 

(5.12) logPie'^'') ^ -XK^hHog{l/{\K\h)) if a = 2, k^O. 



Recalling that it suffices to prove (3.11) for k 7^ and observing, that there 
the parameter k. can be treated like a constant, we sec that log(l/(|/«|/i)) ~ 
log (l/h) , where, because h ^ , we can assume < h < 1 . 
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Hence we can replace (5.12) by 

(5.13) logp(e^'^'*) ~ -XK'^h? log ^ if a = 2, k^Q. 

f b 

Then the limit relation (3.11) (equivalently (3.10)) holds if we scale by 

(5.14) r = (7(/i) = — -^^!— — /i" if 0<q;<2, 

^ ' ^ ^ r(a + 1) sin(a7r/2) 



(5.15) T = a{h) = \h? log ^ if a = 2 . 



Putting ji — A7r/(r(Q; + 1) sin(Q;7r/2) = X/b{a) in (5.9) with b{a) defined in 
(2.17) we obtain from (5.4) the regular scahng law 

(5.16) T = a{h) = fih"^ for < a < 2 
with the restriction (5.6) for fi . As result we have 

Theorem 5.1. Distinguish the cases (i) < a < 2, (ii) a = 2. Define 
the probabilities pk = P{Y = k) in case (i) by (5.4) with restriction (5.6), in 
case (ii) by 

Po = l-2AC(3), Pk^X\k\-^ for k^O 
with restriction < A < 1/(2(^(3)) . Let the scaling relation 

(5.17) T = iih°' in case (i), t — \h°'\og— in case {ii) 

hold and let for fixed t > the index n = t/r run through N towards oo . 
Then the random variable Sn of (3.3) converges in distribution to the random 
variable S{t) whose probability density is given by (1-15) as ga{x,t\Q) . 
Remcirk 5.2. We can use throughout < A < 2 the parameter A and 
then have in (5.1) under the restriction (5.2) a unified representation of the 
transition probabilities. Here, in contrast to the Griinwald-Letnikov random 
walk, the value a = 1 does no longer play a special role. With = X/b{a) we 
have for < a < 2 the regular scaling law r = /ih"' . However, the price to be 
paid for this unified representation is the non-regular scaling t — log(l//i) 
for a = 2 . Another price is that the generating function p{z) in (5.8) is non- 
elementary, requiring considerable efforts in its asymptotic analysis. 
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6 A globally binomial random walk 



The random walk model discussed in Section 4 has the disadvantage that the 
case q; = 1 is excluded and the representation of the transition probabilities 
Pfc for 1 < a < 2 is different from that for < a < 1 . However, for all 
admissible values of a we have the regular scaling law r = ///i". The method 
treated in Section 5 has the advantage of a unified representation of the 
transition probabilities in the whole interval < a < 2 , but the scaling law 
T = iih" holds only for < a < 2 , it breaks down at a = 2 . In this section 
we present a model that in the whole interval < a < 2 admits a unified 
representation of the via binomial coefficients and has there a scaling law 
of the form t = /j.h"' . Moreover, the generating function p{z) is elementary 
for all a E (0, 2] . 

The use of the binomial coefficients in the Griinwald-Letnikov random 
walk has caused singular behaviour for a = 1 . One reason for this sad fact 
is that = for integer j > 2 . We can remove this singular behaviour by 
removing the factor a — 1. 

For < q; < 2 , a 7^ 1 let us define 

(6.1) po = 1 - 2A, pk = (-1)'+'^ d^i"^ ^ ioT k^O. 

Observing that here the singularity at a = 1 is removable, let us for a = 1 
define (via a — )■ 1 in (6.1)) 

(6.2) p. = l-2A. P. = ^^(p^ for MO. 

In (6.1) and (6.2) JZ^^ZPf^ ~ ^ if < A < 1/2 all > . In the special 
case q; = 2 we get 

Po = l-2A, pi=p_i = A, pfc = for |A;| > 2, 



the familiar random walk for approximation of the classical process governed 

du ^ d^u 
dt dx'^ 



by the equation ^ — 



The generating function p{z) = Y^^^^'^Pk^ has in the case a ^ 1 the form 
(6.3) p{z)^l-\{q{z) + q{z-^)} 
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with 

q{z) = ^^{l-z-'){{l-zr-'-l}. 

By passing here to the hmit or directly from (6.2) we get for a = 1 the 
representation 

(6.4) p{z) = 1 - A{(1 - z-^) log(l -z) + {l-z) log(l - z-^)] , p(l) = 1 . 

We have proposed and investigated the particular random walk so generated 
(its transition probabilities given in (6.2)) in [GM99, Section 5]. 

In the special case ct = 2 we find 

(6.5) p{z) = l + X{z-2 + z-^). 



We will now show that for all a e (0, 2] there exists a finite positive number 
c{a) so that, with 

(6.6) II — c{a) A , 

we arrive for e R \ {0} at the small h asymptotics 

(6.7) Pie'^^) = 1 - //(|«|/i)° + o ((|«| /i)") 

which implies (3.11). As in Sections 4 and 5 we can ignore the value n = 
as trivial. 

Referring to [GM99] for detailed treatment of the case o; = 1 , let now be 
7^ K e R and 0<q;<2, a ^ 1 , z — e*'*'*. In view of (6.3) we investigate 
the asymptotics of q{z) + q{z~^) for /i — >■ 0. Prom z~^ — z and 

(1 - a)q{z) = z-'{l - zY - + 1 = e-'''^{\ - ^''^Y - e''"^ + 1 , 

we conclude on 
(6.8) 

^P{z) := (1 - a){q{z) + q{z-^)} = 23? {e-*''^(l - e"*'')"} + 2(1 - cos(«:/i)) , 
and here 

(6.9) ^ {e-'^''{l - e'^^'Y] - 3? ((-m h^) = h)" cos (a7r/2), 
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(6.10) 1 - cos{Kh) l-{\K\hf . 



We distinguish three cases: (i) < a < 1 , (ii) 1 < a < 2 , (iii) a = 2 . 

In cases (i) and (ii) the leading term in the asymptotics of '4){z) turns out to 
be 

i){z) ~2(|fi;|/i)" cos(a7r/2). 

In case (iii) where a = 2 however, this term is matched in order of magnitude 
by (6.10) so that we obtain 

i){z) - 2(|av| hf{-\) + {\k\ hf = hf . 



Collecting results and dividing (6.8) by 1 — a we get (with z — e*'*'*) 

, 2 cos (ail 12) , , , 

,A ^ ^(k/i" ifO<Q;<2,Q;7^1, 

(6.11) X{q{z) + q{z-^)} r-. { 1 - a ' ' 

X{\K\hf if a = 2. 



Hence, in view of (6.3), we obtain (6.7) with (6.6) by putting 
(6.12) c(a) = < 



f 2cos(a7r/2) 

^ ' ' if < a < 2, a ^ 1 , 



1 — a 

1 if a = 2 . 

The scaling coefficient c{a) allows continuous extension to the value a — 1 , 
giving lim c(a) = vr in accordance with [GM99, formula (5.1)1. At a = 2 , 

a— >1 

however, c{a) is discontinuous. In fact 

(6.13) c(2) = 1 7^ 2 = lim c{a) . 

a— ^2 



Let US finally display the transition probabilities with instead of A as pa- 
rameter. 
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For < q; < 2, a^l: 



(6.14) 



a. 



Pk = 



2 cos(Q;7r/2) ' 



0< H< 



2cos(Q;7r/2) V|A;| + 1^ 
cos {an/ 2) 



for ky^O, 



1 — a 

for a = 1 (see [GM99, formula (5.1)]): 



(6.15) 



2/i 

Po = 1 , Pk^ - 

TT TT 



2: 



+ 1) 



for A; 7^ , < // < 7r/2 , 



for a; 
(6.16) 

Po = l-2//, pi^p_i^i_i, pfc = for |/i;|>2, < < 1/2 . 
The discontinuity at o; = 2 has so been transferred to the upper bound for 

We comprise the resuh in 

Theorem 6.1. Take the probabilities Pk = P{Y = k) and the restrictions 
for IX as in formulas (6.I4), (6.15), (6.16), and use the scaling relation 
T = fih'^ . Let for fixed t > the index n — t/r run through N towards 00 . 
Then the random variable Sn of (3.3) converges in distribution to the random 
variable S{t) whose probability density is given by (1.15) as ga{x,t;0) . 



7 The Chechkin-Gonchar random walk 

In this section we adopt to each other considerations of Chechkin and Gon- 
char [ChG] and the framework of our Section 3, restricting attention to the 
parameter range < a < 2 . So doing we exclude the well-known case of the 
classical Gaussian process. We will obtain a random walk, which is discrete 
in time but continuous in space, in more precise words: whose jumping width 
(in the instants tn — nr) can assume any real number, having an everywhere 
positive probability density. We modify our theory of Section 3 by allowing 
the random variable Y to have a strictly monotonic continuous distribution 
function W{x) — P {¥ < x) {x e R), that we furthermore require to be 
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symmetric in the sense W{x) + W(—x) — 1 (x G R). being only interested 
in the symmetric case ^ = of Section 1. We then have in the sum 

n n 

(7.1) ^ (/irj = /i (neN) 

m=l m=l 

a description of a random walk, starting in the point x — . Here h > 
is a scaling width that we let depend on the time-step r > via a strictly 
monotonic scaling relation r = a{h) , with a{h) — as /i — )■ . We expect 
the scaling relation to have the form r = fih°' with the positive coefficient 
jj, to be specified, by having found orientation in Gnedenko's theorem on 
normal attraction (see [GnK], §35). It should be noted, however, that in this 
theorem the scaling constant C appearing there is given with a wrong value 
as has been remarked in [Ba]. 

As previously, we let the Ym be independent identically distributed random 
variables, all having their distribution common with Y . However, we now 
assume Y to have an everywhere positive (not necessarily bounded) proba- 
bility density w — W which is an even function w{x) = w{—x) (x G R) . 
We will use the fact that w is normahzed, w{x) dx — 1 . 

Fixing a value i > and again setting t — tn — nr (equivalent to n — t/r) 
with n G N we want that the random variable Sn converges in distribution to 
the random variable S{t) whose density is given by (1.15). To this purpose we 
introduce a condition on the asymptotic behaviour of the density w, namely 
(7.2) 

w(x) = (6 + e(|x|)) \e{\x\)\<mm [k,E\x\-^} (x G R), 

with positive constants b , K , E and 7 . 

With w{k,) = Jl^^ e^'^^ w{x) dx as characteristic function of the density w{x) 
we observe that the random variable hY has density w{x/h)/h, hence the 
characteristic function w{Kh) , and proceeding in analogy to the general 
method described in Section 3, replacing p{k,, h) = p{e^'^^) in (3.9) by w^nh) , 
we will find a scaling function cr(/i) such that for all k G R, in analogy to 
(3.11), 

(7.3) \ogw{Kh) ^ as h^Q. 

a[h) 

Of course, (7.3) is trivial for k = . Since w like w is an even function 
it suffices to consider (7.3) for (fixed) values ac > . We see that (7.3) is 
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equivalent to 

(7.4) w{Kh) = 1 - a{h) + o{a{h)) as h^O. 

In view of the symmetry and normahzation properties of w {x) and abbrevi- 
ating nh = u we find 

P OO /"OO 

w(u) - 1 = (e'"''' + e-*^^ - 2) w(x) dx = -4 jj^ (sin(i/x/2))^ 
so, using (7.2), 

w{u) = l-2-"+Hi/" / r""^ (sinO / e(x) x"""^ (sin(i/x/2))^ dx . 

Jo Jo 

The first integral can be evaluated in terms of the gamma function. In fact, 
from [GR, (3.823)] we take 



roo 

Jo 



r(-n) cos(nA/2) 



2i-a 22-" r{a + 1) sin(Q;7r/2) 

The latter equality follows by the reflection formula for the gamma function. 

We estimate the second integral via decomposition /q°° . . . = Jq . . . + . . . , 

taking 7] = z/-(2"+7)/{2a+27)^ using |sin^| < min{^, 1} for ,^ > and the 
condition on e(|a;|) of (7.2). By careful calculation we find that it behaves 
asymptotically as o (i/") = |«;|" o {h^) . Combining these results and recalling 
that w is an even function, we obtain 

(7.5) wUh) ^ I - \kV' — --/i" + kro(/i") (/i^O) 

^ ' ^ ' ' ' r(a + 1) sin(a7r/2) \ \ \ ) \ ) 

as valid for all k G R . In view of (7.4) and the theory developed in Section 
3 we thus arrive at the scaling relation 

(7.6) T^a{h)=iih", with ji 



r{a + 1) sin(a7r/2) ' 



Now we are in the position to formulate 

Theorem 7.1. Let < a < 2 and assume the random variable Y to have 
a probability density w of the form (7.2). Let the scaling relation (7.6) hold 
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and let for fixed t > the index n — t/r run through N towards oo . Then the 
random variable Sn of (7.1) converges in distribution to the random variable 
S{t) whose probability density is given by (1.15). 

Remark 7.2. According to the well-known asymptotic expansions of the 
function Pa (a;; 0) = ga{x, 1;0) (see [F52], [F71], [Zo]) we have 

^7 7^ ^ ^ r(Qi + 1) sin(Qi7r/2) 

TT 

if we take w{x) — Pa{x; , hence in this case t — 1 , /j, = 1 and h = r^/" = 
^-i/q ^ j£ require in (7.1) the Ym to have this special density, then Sn for 
all t > has the same probability distribution as S{t) whose characteristic 
function is exp(— t \k\"') . We can here obtain the scahng relation also via the 
convolution theorem. 

Remcirk 7.3. For actual simulation a random variable Y having the re- 
quired properties is particularly useful if its distribution function W{x) = 
f!^^ w{C,) dC, is easily invertible. We can then generate a realization of F by a 
standard Monte Carlo method (see [HH]). Generate a random number y uni- 
formly distributed in the interval [0, 1) . Then solve the equation y — W{x) 
for X and take realization of Y. Chechkin and Gonchar in [ChG] have 

proposed to use 



(7.8) W{x) 



i (1 + |a;|")"^ for x <0, 

1 - i (l + x°)~^ for x>0, 
a function easily invertible. The density 

(7.9) w{x) = W'(x) 





\x[ 


a— 


1 


2 (1 + 


\x\ 


aj2 



has the property (7.2) with b — a/2 , ^ — a , hence we get 



(7.10) /X 



TT 



2r(Q;) sin(Q;7r/2) ' 



The density (7.9) is unbounded at the origin if < a < 1 . To avoid this we 
propose 



(7.11) w{x) = -{l + \x 
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which again satisfies the asymptotic condition (7.2). Then 



(7.12) 



W{x) = 



for X < , 



I - ^ (l + x)"" for x>0 
is also easily invertible, and (7.10) remains valid. 

Remark 7.4. Among the symmetric densities Pa{x] 0) , only the Cauchy 
density w{x) = pi{x;0) = (I/tt) (1 + x'^)~^ offers easy invertibility of the 
corresponding distribution function, namely of the function W{x) = 1/2 + 
(I/tt) arctanx . Via random numbers Um uniformly distributed in [0,1) we 
can get reahzations of the in Tta,n{7r{yjn — 1/2)) (here t — Ih^ — h) and 
so obtain in Sn a snapshot at instant tn — nr of a true Cauchy process. 



8 Conclusions 

Anomalous diffusion processes have in recent years gained revived interest 
among physicists, and methods of fractional calculus have shown their useful- 
ness for purposes of modelling. In the space-fractional case one is naturally 
led to a generalization of the classical diffusion equation with respect to the 
second-order spatial operator. One arrives in a natural way at the processes 
of Levy-Feller type in which stable probability distributions play the essential 
role. Also among physicists and mathematicians who have found it reward- 
ing to work in theory of finance, such processes are becoming more and more 
popular (see e.g. [M], [BoP], [MS]). So, it is no wonder that also in pure 
mathematics such types of processes are now investigated in great generality 
and analytical sophistication (see e.g. [J], [Be], [S], [Za]). From the more 
practical point of view discrete models are esteemed. They not only show 
that very different microscopic behaviour of particles can result in the same 
macroscopic behaviour but offer also possible visualizations of what is hap- 
pening in such processes. Furthermore such discrete models can be used for 
simulation purposes, be it for simulation of particle paths via Monte Carlo 
methods (the microscopic view) or via solution of the underlying Cauchy 
problem for a pseudo-differential equation (the macroscopic view). And, last 
but not least, such models are fascinating as seen from the mathematical 
standpoint (or, more specifically, from the position of probability theory). 
In our present investigation we first have given a survey on and drawn mo- 
tivations from basic theory of fractional calculus and Levy-Feller diffusion 
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processes. Then we have obtained and rigorously analyzed (with respect to 
their convergence in distribution for passing to the limit of infinitely fine 
discretization) three models of random walk occurring on a regular spatio- 
temporal grid. The first model is devised from the Griinwald-Letnikov dis- 
cretization of the two Weyl operators, the composition of which gives the 
inverse of the Riesz potential operator. The second model is an adaptation 
of ideas of Gillis and Weiss [GiW] to our framework. We have provided it with 
a new motivation, namely as obtainable from straightforward discretization 
of the hypersingular integral representation of the spatial pseudo-differential 
operator. The third model's intention is to overcome peculiar deficiencies of 
the first two models. It is a modification and improvement of the first model, 
and again properties of the binomial coefficients are used. 

Finally, to offer also a highly efficient method for numerical simulation, we 
have mutually adapted our theoretical frame to ideas of Chechkin and Gon- 
char [ChG]. We so obtain a random walk still proceeding in equidistant in- 
stants of time but allowing spatial jumps of arbitrary length in positive or 
negative direction. 



Appendix A: Asymptotics of an integral 

Abbreviating nh — u in z — e^'^^ in the right hand side of (5.10), and keeping 
in mind < a < 2 , elementary calculation yields the equation 



(Al) T{a + l)^^{z) = r 

Jo 



(1 - e-") |1 - e-"e' 



iv\2 



which we will treat asymptotically for < v ^ 0+ by the Laplace method for 
integrals (see [dB]), using the fact that the lower bound w = is the critical 
one (the integrand tending to oo as — >■ 0). We have 1 — cosi/ = u'^/2 + 0{u'^) 
and 

|1 _ e-^e^'^l' = (1 - e-")2 + 2e-"(l - cosu) = (1 - e-")^ + i/^e"" + 0{u^), 
uniformly in < < oo, hence 

z/2 foo w"e-"(l e-") 



du . 
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Because this integral diverges for v — we can simplify the integrand (for 
small u) which, for small u, gives the essential contribution: 1 + e~" ~ 2 , 
1 — e~" ~ u , e~" ~ 1 . We obtain 



r(a + 1) 3?7(z) ~ z/^ / 
Jo 

and, by substituting u — uw, 

iA.2) Via + l)^-^{z) = i/" / w"-^^^ dw := . 

^0 + 1 

In the investigation of the integral 

(A3) pM = / -"-^d^. 

we distinguish the cases (i) < a < 2, (ii) o; = 2 . In the case (i) simply 

p{^) ~^ / ~^ for i/ — )■ 

Jo w-^ 

and with ^ = a — 1, hence — 1 < ^ < 1, we have to determine the value of 

fOO rpP 

Observing that q{—l3) = q{(5) (substitute ^ = 1/x) we do this < ^ < 1 . 
Complementation by (integrate along the upper edge of the negative real 
semi-axis) 



-oo + 1 Jo X'^ + 1 

gives, via the residue theorem, 

(X r+oo rpP 

1 + e^^'^) g(/3) = / ^— - dx^ 1:1^^1: 
' J —00 X -\- \ 

So 

m 



2cos(^7r/2) 2sin(Q;7r/2) ' 
and hence 

(^•^) "M^iita^ if 0<a<2 (.^0+). 
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In case (ii) the integral diverges for z/ = 0, so we must proceed in another 
way. Inserting a = 2 in (A. 3) and differentiating we obtain for i/ > 

Jq w^-\-l Jo \ + 1/2 

and then by integration 

{A.5) p{v) ~ - log i/ = log - {v ^ 0+). 



Now we can collect results. From (A.l) - (A. 5), using v = nh which because 
of symmetry we can replace by \K\h (admitting also negative values of k) we 
deduce 
(A.6) 

^7(-g)- ^p, ^-,7- ( if 0<q;<2, «^0, as /i ^ , 

2T[a + 1) sm [aTi/2) 



(A7) 3f?7(^) ~ /«2/iMog ^ if a = 2, , as h^Q. 

\K\h 
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